Ideal evolution of MHD turbulence when imposing Taylor-Green symmetries 



M.E. Brachet 1 , M. D. Bustamante 2 , G. Krstulovic 3 , P.D. Mininni 4 ' 5 , A. Pouquet 4 and D. Rosenberg 4 
1 Laboratoire de Physique Statistique de I'Ecole Normale Superieure, 
associe au CNRS et aux Universites ParisVI et VII, 24 Rue Lhomond, 75231 Paris, France. 
2 School of Mathematical Sciences, University College Dublin, Belfield, Dublin 4, Ireland. 
3 Laboratoire Lagrange UMR7293, Universite de Nice Sophia- Antipolis, CNRS, 
Observatoire de la Cote d'Azur, B.P. 4229, 06304 Nice Cedex 4, France 
4 Computational and Information Systems Laboratory, 
NCAR, P.O. Box 3000, Boulder CO 80307, USA. 
5 Departamento de Fisica, Facultad de Ciencias Exactas y Naturales, 
Universidad de Buenos Aires and IF IB A, CONICET, 
Ciudad Universitaria, 1428 Buenos Aires, Argentina. 

We investigate the ideal and incompressible magnetohydrodynamic (MHD) equations in three 
space dimensions for the development of potentially singular structures. The methodology consists 
in implementing the four-fold symmetries of the Taylor-Green vortex generalized to MHD, leading 
to substantial computer time and memory savings at a given resolution; we also use a re-gridding 
method that allows for lower-resolution runs at early times, with no loss of spectral accuracy. One 
magnetic configuration is examined at an equivalent resolution of 6144 3 points, and three different 
configurations on grids of 4096 3 points. At the highest resolution, two different current and vorticity 
sheet systems are found to collide, producing two successive accelerations in the development of small 
scales. At the latest time, a convergence of magnetic field lines to the location of maximum current 
is probably leading locally to a strong bending and directional variability of such lines. A novel 
analytical method, based on sharp analysis inequalities, is used to assess the validity of the finite- 
time singularity scenario. This method allows one to rule out spurious singularities by evaluating 
the rate at which the logarithmic decrement of the analyticity-strip method goes to zero. The result 
is that the finite-time singularity scenario cannot be ruled out, and the singularity time could be 
somewhere between t = 2.33 and t — 2.70. More robust conclusions will require higher resolution 
runs and grid-point interpolation measurements of maximum current and vorticity. 

PACS numbers: 47.10.A, 47.65-d,47.15.ki,47.11.Kb 



I. INTRODUCTION 

The class of problems addressing the formation of sin- 
gularities and the existence and structure of solutions 
of nonlinear partial differential equations for all times 
forms an important branch of mathematics, with wide 
application in numerous fields: engineering, astro- and 
geophysics, laboratory studies of superfluids, and in me- 
teorological research on extreme events such as torna- 
does and hurricanes. The presence or absence of either 
dissipation-viscosity, of magnetic resistivity in magneto- 
hydrodynamics (MHD), or of dispersion, plays an essen- 
tial role as well. The significance of such questions is 
recognized for example by the Clay Institute Millennium 
Prize for a proof of existence and smoothness of finite en- 
ergy solutions of the Navier-Stokes equations, and by the 
numerous studies devoted to them: How fast do (poten- 
tially) singular structures form? What is their temporal 
evolution and geometry? What role do their interactions 
play and how might they lead to a modification of trans- 
port properties within complex flows, including in the 
presence of magnetic fields? Progress on such problems 
will most likely come from a combination of mathematics, 
laboratory experiments and direct numerical simulations 
(DNS), in the latter case in particular using computer 
codes with high accuracy and performing studies at the 
highest possible resolutions. 



There is a large body of analytical and numerical work 
on singularities in fluids. As theoretical estimates are 
not necessarily sharp, numerical data are invaluable in 
assessing potential singularities, as discussed, e.g., in pQ. 
Unfortunately, in the case of the numerics, with regard to 
existence the answer vacillates between "yes" and "no" 
[2]. Singularities occur in simplified models, as derived 
in [31 3] assuming an isotropic pressure Hessian. These 
models have been generalized to MHD in the vicinity of 
magnetic null points [5 and lead as well to a singularity, 
but the question remains open in the general (and most 
physically relevant) case. 

One of the most useful criteria in the search for a singu- 
larity comes from the Beale-Kato-Majda theorem (BKM 
hereafter) [6] which states for incompressible ideal fluids 
that, if the flow presents a finite-time singularity at T*, 
then 

/ ||w(.,t)|| 0o £ft = C», (1) 

Jo 

where we have used the usual notation for the or 
supremum norm, u> = V x v being the vorticity and v 
the velocity. If a power-law divergence of vorticity at T* 
is assumed, ^(.,1)^ = C\T* - t\~ , for t -> T«, with 
/3 > and C a constant, then the BKM theorem can be 
re-expressed as: "The flow has a finite-time singularity at 
t = T* if and only if (3 > 1." As stressed in [7j, enstrophy 
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production Dtt/Dt, with fi = J a; 2 (x)<i 3 x (the £2 norm), 
should also be monitored to detect singularities, and care 
must be taken in assessing C, (3 and when fitting the 
data stemming from DNS. 

Furthermore, the dynamics of the vorticity (and cur- 
rent in MHD flows) should be monitored not only for 
their C2 and norms, but also for changes in the di- 
rection of their field lines (or "swing" [8 ). It was found 
in neutral flows that the rapid growth of ||o;|| can be 
countered by the straightening of vortex lines (see [8THU] 
in the MHD case). Moreover, the study of the evolu- 
tion of the curvature and torsion of vortex (or current) 
lines yields interesting insights into the dynamics of ideal 
flows pp. The rich variety in the observed behavior has 
resulted in a plethora of initial conditions examined in 
previous numerical studies. 

Among the 3D flows that have been considered for their 
potential singular behavior for ideal (non-dissipative) flu- 
ids are the Taylor-Green flow (TG hereafter) [11], the 
Kida-Pelz flow (KP) [T2J[T3], an d two anti-parallel vor- 
tices P3JE3, all displaying symmetries that can be imple- 
mented numerically (see also IB; ). These flows have been 
studied by several teams, with a recent revival IT7H20] 
(see, e.g., [2] for a brief introduction to the literature). 

In MHD when coupling to a magnetic field, the theo- 
rem equivalent to BKM involves the sum of the maxima 
of vorticity u> and current density |21j. Ideal MHD in two 
space dimensions has been studied in the past (see [2"2T 
and more recently using high-resolution runs |27|). 
but in the 3D general case, ideal runs are scarce except 
for the pioneering work using symmetric configurations 
of linked flux tubes with zero initial velocity |28| , or with 
adaptive mesh refinement (AMR) using finite differences 
if]. 

One can also use the TG flow and generalize it to MHD 
(hereafter, TG-MHD flows), as done in [32]. One of the 
TG-MHD flows studied for its possible singular behav- 
ior in [32 displays a feature not observed at the time in 
the fluid case: after an initial phase of thinning of the 
current and vortex sheets, the flow outside the structure 
pushes together two current sheets with widely different 
directions of the magnetic field embedded in them, lead- 
ing to a rotational quasi-discontinuity with a substantial 
acceleration in the development of small scales. Once dis- 
sipation is restored, this small-scale activity is diagnosed 
as intermittent reconnection 33 . Rotational and tangen- 
tial discontinuities, identified as intermittent structures, 
have been observed in the Solar Wind using a variety 
of in-situ acquired data [51]; they have also been iden- 
tified at the edge of Reverse Field Pinch plasma devices 
(see, e.g., [3S] for review). Using the Cluster ensemble 
of four satellites, all four spacecrafts indicate at times a 
directional (either rotational or tangential) discontinuity, 
including with a small normal component of the magnetic 
field. Such rotational discontinuities can stem from non- 
linear steepening or from reconnection of magnetic field 
lines [37]. Their modeling leads to statistical properties 
akin to that of so-called nano-flares observed in the so- 



lar corona [3J5], and they provide tantalizing hints that 
singularities may exist in MHD. 

Only by performing substantially higher-resolution 
and high accuracy runs that high-performance comput- 
ing resources can allow, shall we be able to explore sev- 
eral configurations leading to possible singular behavior 
in MHD. It is in this context that we propose to search 
in this work for singularities in MHD with different con- 
figurations and using the highest known resolutions (and 
hence, scale separation between the size of the box and 
the size of the mesh); thus, a run is performed on an 
equivalent grid of 6144 3 points in one case, following on 
the work done in [32] on grids of 2048 3 points. 

II. NUMERICAL PROCEDURE 

In this section we describe briefly the codes and meth- 
ods used, and give details of the numerical simulations. 
We present first the MHD equations, and then intro- 
duce the initial conditions. Then, we explain how the 
code is parallelized for the simulations at the largest res- 
olutions. The choice of de-aliasing method is crucial to 
conserve the total energy and other quadratic invariants 
with good accuracy, and details concerning our method- 
ology are given and compared with other choices for de- 
aliasing. Then, the procedure followed to increase spatial 
resolution as structures become thinner is explained. Fi- 
nally, we comment on the effect that imposing the four- 
fold symmetries of the Taylor-Green vortex generalized 
to MHD might have. 

A. Equations, initial conditions and the TYGRS 
code 

The MHD equations for an incompressible and ideal 
fluid with v and b respectively the velocity and magnetic 
field read: 

^+v-Vv = -— VP+jxb, (2) 
^ = V x (v x b) ; (3) 

Po = 1 is the (uniform) density and b is the Alfven veloc- 
ity, V is the pressure, V-v = 0,V-b = 0, and there are 
no dissipative or forcing terms; finally, j = V x b is the 
current density. The total (kinetic plus magnetic) energy 
Et, the cross helicity He and the magnetic helicity Hm, 
defined as 

E T = E V +E M = (v 2 + b 2 )/2 (4) 
H c = {vb) , H M = {A-b) (5) 

with A the magnetic potential (b = V x A), are all 
conserved by the nonlinear interactions (38] - 

In practice, a pseudo-spectral code solves these equa- 
tions in Fourier space, truncated up to some maximum 
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wavenumber. The truncated MHD equations for the 
Fourier modes Uk and bk, with k € [k m i n ,k max ) can be 
written easily, the Fourier modes satisfying u k = 0, b k = 
if |k| > /c max or if |k| < fc m i n . For a computational box 
of length 2-7T, we have /c m in = 1, and with a de-aliasing 
using the 2/3-rule, fc max = A/3, where N is the num- 
ber of modes per dimension (we assume a box with unit 
aspect ratio). Other de- aliasing methods can be used 
succ essfu lly [P71 \TE\ and are discussed briefly below (see 
Sec. II D[ ). It is important to note here that de-aliasing 
is crucial in pseudo-spectral simulations to remove spu- 
rious growth of modes with large wavenumbers, and to 
conserve the total energy and other quadratic invariants. 
Indeed, a pseudo-spectral code which is fully dealiased is 
equivalent to a Galerkin truncation, and thus preserves 
all quadratic invariants in the system to round-off error. 

The equations are solved starting from initial condi- 
tions for the velocity and the magnetic field. If the ini- 
tial conditions have symmetries that are preserved by 
the equations, then the symmetries can be used to save 
memory and computing time. As already mentioned, in 
hydrodynamics (b = 0) one of the simplest velocity fields 
satisfying the symmetries of the equations is the TG flow 
(note that the z component, initially equal to zero, will 
grow with time) [TTJ \M SD]: 

u(x, y, z) = uq [(sin x cos y cos z)e x — (cos x sin y cos z)e y ] 

(6) 

It is interesting to point out that the TG flow in a pe- 
riodic box shares similarities with the von Karman flow 
between two counter-rotating disks as used in several lab- 
oratory experiments, including those with liquid metals 
such as sodium or gallium, to study the generation of 
magnetic fields. 

To generalize the TG flow to MHD, we use the veloc- 
ity as prescribed by Eq. ([6]), and we will consider three 
possible choices for the initial magnetic field b with the 
same overall symmetries [3H [33J EI] • We refer to these 
three flows as the insulating (I) defined by 



cos x sin y sin z 
b 1 = &n I sin x cos y sin z 

-2 sin x sin y cos z j 

the alternative insulating flow (A): 

cos 2x sin 2y sin2z 
6q I — sin 2x cos 2y sin 2z 


and the conducting flow (C): 

sin 2x cos 2y cos 2z 
cos 2x sin 2y cos 2z 
—2 cos 2x cos 2y sin 2z, 



(7) 



(8) 



(9) 



Note that the I, A and C flows, with almost identical in- 
variants, have nevertheless three different developed en- 
ergy spectra in the non-ideal case at the maximum of 



dissipation [33J , displaying a lack of universality in MHD 
turbulence in the absence of an imposed magnetic field. 

For all three confi gura tions, Ey = Em = 0.125 when 
u = 1 and b = \/l/3, 1, y/2/3, respectively; for the 
helicities, Hm = and He ~ because of the imposed 
symmetries; note however that there can be strong local 
correlations corresponding to local alignment of u and 
b, as can be shown both analytically and numerically 
[43] . In the I case, the current j l is everywhere paral- 
lel to the walls of the so-called impermeable box [0,7r] 3 
which thus appears to be insulating. For the C case, j c 
in the [0, 7r] 3 box is perpendicular to the walls, which are 
therefore conducting. In this configuration, He is non- 
zero but small (less than 4% at its maximum over time, 
in a dimensionless measure relative to the total energy). 
Finally, b a is an alternative insulating MHD vortex. 

The code, TYGRS (TaYlor-GReen Symmetric; see be- 
low), enforces the symmetries of the TG vortex in 3D 
hydrodynamics, and of the TG-MHD vortices in 3D 
MHD within the periodic cube of length 2ir. These 
symmetries include: mirror symmetries about the planes 
x = 0&7r, y = 0& 7r and z = & tt together 
with x — tt, y — 0, y — tt, z — 0, and z = tt 
(e.g., in the x-direction: v x {— x,y, z) — —v x (x,y,z) and 
v x (TT — x,y,z) = —v x (-TT + x,y,z)), rotational symme- 
tries of angle rnr about the axes (x,y,z) = (§ ,y, f) and 
(x, £ , §), and rotational symmetries of angle mr/2 about 
the axis (f , ? ,z), for n e Z. Because of these symme- 
tries, the Fourier-transformed fields are non-zero only for 
wavenumbers (k x , k y , k z ) with jointly even or jointly odd 
components. 

Thus, TYGRS computations at a given scale separa- 
tion (defined as the ratio k ma x/kmin, which is propor- 
tional to the Reynolds number in the dissipative case), 
or at a given equivalent resolution, are performed on lin- 
ear grids that are one-fourth the size of those for a general 
code, by exploiting symmetries of the TG vortex: one ob- 
tains the flow in the full periodic box of size [0, 2tt) 3 by ap- 
plying these symmetries to the impermeable box [0,7r] 3 . 
The nonlinear terms and their temporal derivatives are 
computed from the even-odd decomposition of the fields 
in the fundamental box [0, tt/2] 3 . Note that TYGRS per- 
forms a DNS, since no modeling of small scales is done. 
For time integration, an explicit 2™ -order Runge-Kutta 
scheme is used. Because the time integration truncation 
error at the proposed resolutions may exceed the single 
floating point precision, we use double precision for the 
computations. 

No uniform external field B will be imposed in our 
simulations. Such an external field is known to slow- 
down small-scale development, and may quench the de- 
velopment of singularities [22j |44] because of the semi- 
dispersive nature of the problem, with Alfven waves prop- 
agating in opposite directions along B . This slowing- 
down of nonlinear dynamics due to waves has been mod- 
eled phenomenologically in several ways, starting with 
Iroshnikov and Kraichnan in the mid sixties with a /c~ 3 / 2 
total isotropic energy spectrum, as opposed to the clas- 
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FIG. 1: Timings of the TYGRS code on Jaguar up to ~ 10 4 
cores and up to an actual number of computational degrees 
of freedom of N% d = 2048 3 points. 



sical Kolmogorov spectrum for fluid turbulence. It can 
be evaluated analytically using weak turbulence theory 
for large Bo [U3 US] , leading to a steeper and anisotropic 
spectrum ~ kj 2 , with k± referring to the direction per- 
pendicular to B . 



B. The role of symmetries 

In [35] , simulations with and without imposed symme- 
tries with Taylor-Green initial conditions were compared. 
No differences were observed except at the lowest mode 
and at an energetic level close to round-off error. Also, vi- 
sualization analyses showed that the physical structures 
that are present in the flow appear identical between the 
runs with and without imposed symmetries (see, e.g., 
[4"T]). Of course, at late times instabilities develop in- 
duced by noise due to accumulated errors because, e.g., 
of insufficient numerical accuracy [42] . These errors can 
break the symmetries in the computation of the flow and 
field, when one does not impose the symmetries of the ini- 
tial conditions. In that case, magnetic and cross-hclicity 
grow and may lead the flow to another final state. How- 
ever, this bifurcation in behavior happens at a signifi- 
cantly later time than the times considered in the present 
study. 



C. Implementation of the hybrid scheme for the 
TYGRS code 

Pseudo-spectral codes are known to be optimal on 
periodic domains |47j . However, they require global 
spectral transforms, and thus are hard to implement 



in distributed memory environments, a crucial limita- 
tion until one-dimensional domain decomposition tech- 
niques (DDT) arose, that allowed computation of serial 
Fast Fourier Transforms (FFTs) in different directions in 
space (local in memory) after performing transpositions. 
However, distributed parallelization using the Message 
Passing Interface (MPI) in pseudo-spectral codes is lim- 
ited in the number of processors that can be used, unless 
more transpositions are done per FFT (thus increasing 
communication). The hybrid (MPI-OpenMP) scheme we 
have implemented for a general code builds upon a one- 
dimensional (slab-based) domain decomposition that is 
effective for parallel scaling using MPI alone [H]. In 
the scheme, each MPI task creates multiple threads us- 
ing OpenMP. This method has been extended in TYGRS 
to the sine (cosine) with even (odd) wavenumber FFTs 
needed to implement the symmetries of TG flows, using 
loop-level OpenMP directives and multi-threaded FFTs. 

The resulting quasi-linear scaling up to ~ 10 4 cores 
for TYGRS, particularly at high resolution, is displayed 
in Fig. [I] The hybrid scheme implemented in TYGRS 
was derived from the method developed |48j for a similar 
pseudo-spectral code-Geophysical High-Order Suite for 
Turbulence (GHOST)-in which no symmetries are en- 
forced, that now shows linear scaling up to more than 
98,000 processors on grids of up to 8192 3 points. 

We note that the hybrid scheme used here is not the 
only way in which to decompose the pseudo-spectral grid. 
An alternative is to retain a pure MPI model [49] in 
which the domain decomposition takes the form of "pen- 
cils" and yields a two-dimensional domain decomposition 
among MPI tasks, where OpenMP is not required. This 
technique is also found to scale well to large core counts, 
although large fluctuations in performance are observed 
even within a given processor-domain mapping. The hy- 
brid method offers a two-level parallelization that may be 
more effective in mapping the domain to the hierarchical 
architectures that are now emerging, and better suited 
for environments with multiple cores per socket. The hy- 
brid scheme may also aid in MPI memory problems, in 
that fewer MPI tasks require less buffer memory. This is 
related to the fact that, by reducing the number of MPI 
processes using threads, we reduce not only the number 
of MPI calls, but also the amount of data that must be 
communicated, and hence the size of the MPI buffers re- 
quired. Finally, this also allows us to use parallel MPI 
I/O in environments with tens of thousands of cores, as 
the number of MPI tasks is only a fraction of the total 
number of cores used. 



D. Choice of truncation at high wavenumber and 
the issue of accuracy 

As explained before, one issue to resolve is how best 
to perform the removal of spurious modes with high 
wavenumber, either via a de-aliasing technique using 
the standard 2/3-rule whereby modes are truncated at 
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FIG. 2: Energy spectrum for a hydrodynamic run on a grid 
of 3072 3 points, together with a N=768 3 truncation of that 
dataset using the bootstrap re-gridding down-sizing tool. The 
two spectra are identical (to within machine round-off) up to 
the truncated grid's maximum wavelength, k max = 256. 



2k max /3 where kmax — N/2 is the maximum wavenum- 
ber of the computation on a cubic grid with N points 
on the side, or by multiplying the r.h.s. of the evolution 
equations with a high-order exponential smoothing func- 
tion p(k) = exp[-mi(2fc/iV) m2 ], as proposed in [71 H71- 
I19j with 77i i = ?TJ2 = 36. Using the latter method, more 
Fourier modes are retained in the computation, leading to 
an enhanced scale separation with which smaller scales 
can be reached for a given grid in an ideal flow, and 
thus the computations can in principle be performed for 
a longer time. 

However, when using the second method the exact en- 
ergy conservation in the computations is lost, as can be 
observed in Table HI Also note that the BKM criterion 
given m Eq. for a singularity to occur is based on the 
supremum norm, which is more sensitive to global numer- 
ical accuracy (truncation) [501 15 1 J and numerical preci- 
sion than the £2 measures. Furthermore, it is straight- 
forward to check that exponential smoothing spoils the 
Galilean invariance v(x, t) — > v(x + \Jt, t) — U in the hy- 
drodynamic case. Because of these drawbacks, the 2/3 
de-aliasing rule is used in the following, either in the form 
of the spherical rule (truncation for |k| > N/3) or cu- 
bic rule (truncation for \k x \ > N/3 or |k a | > N/3 or 
|k z | > N/3), as discussed below. 

E. The concept of bootstrap re-gridding 

Besides the constraints given by time stepping errors, 
from previous experience we know that to preserve ac- 
curacy in the computation of spatial derivatives we also 
need to use double precision arithmetic for a grid size at 
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TABLE I: Time evolution of the relative error on energy con- 
servation AE/E for the b = (hydrodynamic) Taylor Green 
initial data at resolution 512 3 , for two types of spectral trun- 
cation. The "Exponential smoothing" method is described in 
the text, while the "2/3-cubic" represents the 2/3 de-aliasing 
rule using cubic truncation of Fourier space. 

or above 4096 3 points. On the other hand, we also know 
that the smallest grid size is only reached slowly (expo- 
nentially in time as long as singularities do not develop) . 
So we propose the following question: Do we need to 
compute from t = to the final time at the maximum 
resolution N that is eventually going to be needed? In- 
deed, at a given linear resolution N± < N, one can com- 
pute until f jv x with sufficient accuracy, as measured for 
example by the logarithmic decrement technique (see be- 
low) . Then, one can restart the run at and compute 
until Tjv 2 with a grid of size N2, with, say, iV2 = 2Ni 
grid points (not necessarily a factor of 2 of course), and 
this process can be rc-itcratcd (m times altogether) until 
we reach the desired resolution N = 2 m Ni, so that only 
the last fraction of the run is done on the largest grid at 
the highest computational cost (in terms of both memory 
and CPU). 

The implementation of the procedure described above 
requires some care when restarts are performed, from the 
point of view of code development because of parallcliza- 
tion of FFTs on grids of different sizes, as well as care- 
ful checking for accuracy for all norms, e.g., £2 but also 
£oo, as needed for singularity tests. However, this "boot- 
strap re-gridding" scheme allows one to save a significant 
fraction of compute time when carrying out the time in- 
tegration at the highest resolution. In the simulations 
presented here, one can estimate a total cost of 1/3 com- 
pared to the full resolution run starting at t = 0. 

It is worth pointing out that the re-gridding scheme 
can also be used to study the dissipative case if one 
chooses to start the run with the last reliable time of 
the ideal run. For forced runs, the extension of the 
methodology is straightforward. But while it may not 
bring about large savings, it might also be useful in cases 
when the flow displays strong signs of intermittent bursts 
followed by long quiescent periods, as for example in the 
case of the stable (nocturnal) planetary boundary layer 
[52] , with the turbulence being related to the presence of 
jets at low altitude. 

In practice, the re-gridding scheme takes a restart 
dataset in physical space, converts to wave-space, and 
then either truncates to reduce resolution (down-sizing, 
useful when performing comparisons with large eddy sim- 
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ulation runs), or else pads (with zeroes) in wave-space to 
increase the spatial resolution. The final step requires 
an inverse multidimensional transform at the new spec- 
tral resolution in order to convert back to physical space 
at the new resolution, so that the data can be used to 
"restart" at the next resolution. The end result of the 
equivalent down-sizing operation is illustrated in Fig. [2j 



III. THE I CONFIGURATION AT HIGH 
RESOLUTION 

A. Implementation of bootstrapping up to an 
equivalent grid of 6144 points for ideal MHD 

The bootstrapping procedure just described can in 
principle introduce errors in the computational procedure 
that breaks the spectral accuracy of the code; hence, we 
show now that this is not the case, provided one is careful 
enough in choosing the time at which the grid resolution 
is increased. In Fig. [3] (top) is given the normalized to- 
tal energy difference (i.e., with respect to initial energy) 
as a function of time, with most of the error occurring 
at early times since the time-step is adapted to the grid 
spacing, which is larger earlier in the computation; the 
different colors (line types) indicate different grid resolu- 
tions. The energy difference remains lower than 10~ 9 at 
all times but shows a rapid increase at the latest times, 
indicative of a build-up of errors. When examining the 
total energy spectra for different times, computed on dif- 
ferent grid resolutions, one can observe a smooth transi- 
tion from one grid to the next (not shown). It is impor- 
tant to note that the re-gridding is performed when the 
energy spectrum at the largest wavenumber in the sim- 
ulation with the grid Ni reaches the machine round-off 
level, with a cut-off conservatively chosen to be 10~ 30 in 
order to preserve a high level of accuracy throughout the 
run. 

Apart from following the numerical conservation of the 
invariants of Eqs. ^ and pi, with special focus on the 
total energy, one diagnostic has been traditionally to 
monitor the logarithmic decrement S, when fitting the 
Fourier spectrum as 



E x (k,t) = cx{t)k- nx ^e-^ Sx ^ 



(10) 



where X stands for either the kinetic (X=V), magnetic 
(X=M) or total (X=T) energy, or the energies of the 
Elsasser variables E± for the fields z± = v ± b. As long 
as 5x =/= the fields remain regular, and when Sx be- 
comes comparable to the mesh the computation of the 
behavior of the partial differential equations ([2| and ([3| 
stops, since at later times one enters the regime of statis- 
tical equilibrium. The logarithmic decrement Sx refers 
to the width of the analyticity strip in the complex plane: 
as long as the complex singularities do not reach the real 
axis, the computation remains regular |44l 153] . Figure 
[3] also gives the temporal evolution of the logarithmic 
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FIG. 3: Top: Normalized energy difference [E T (t) - E ]/E , 
with Eo the initial total energy, showing the error growth 
with time. Note that most runs at a given resolution have 
been pursued for times longer than the time at which re- 
gridding was performed. Error growth is clearly slowed down 
by accuracy (with smaller grid resolution and smaller time- 
step at later times). Different colors (colors online) are used 
for different grid sizes, Ni points per direction with Ni taking 
values: 1536 (blue, solid); 3072 (green, dashed); 4096 (red, 
crosses); and 6144 (black; dash-doted). Middle and bottom: 
Temporal evolution of the logarithmic decrement S and spec- 
tral index n (bottom; see Eq. 10 1 for the total energy spectrum 
(wavenumber fit interval: [10, 1000]). The horizontal line of 
crosses indicate the grid resolution limit 4/Ni for a given com- 
putation on a given grid Gi at an equivalent resolution of Ni. 
The color and line types are the same as in the top figure. 



decrement St (middle) and of the spectral index (bot- 
tom) for the total energy spectrum; grid resolution is 
indicated by the horizontal line of crosses. The fit to 
the spectrum (see Eq. ( 10 ) above) is done in the Fourier 
interval [10,1000]. The acceleration in the decrease of 
the logarithmic decrement found in [32] is confirmed by 
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FIG. 4: (color online) a) Fits (dark or blue lines) from k = 3 
to k = fc max using Eq. ( | 10| ) to the total energy spectra (light 
or red points), in lin-log scale, as a function of wavenumber 
and for different times: t = 1.975, 2.201, 2.425 and t = 2.651. 
Note the good quality of the fit at early times, and the poor 
quality at t — 2.425. b) Same plot in log-log scale. 



the present computation; it is accompanied by a sharp 
increase in the inertial index tit, with both changes oc- 
curring simultaneously at t 2.5. 



However, when comparing the fit using Eq. ( 10 1 to the 



actual spectrum in the simulations, one sees that errors 
are introduced as the spectrum is not always well rep- 
resented by Eq. (10). This is associated with the fact 



that the simple form ( 10 ) needs to be true only in the 
k — y oo asymptotic. We now examine this point fur- 
ther. In simple flows such as the ID Burgers solution 
corresponding to sin(x) initial data, or the purely hy- 
drodynamic Taylor-Green vortex (see |20|). the energy 
spectrum of the flow can be globally well fitted with the 
simple form ( |10[ ), but this is not always the case. For 
instance, in the Kida-Pelz flow, oscillations were found 
and attributed to interferences of complex singularities, 
see [S3]. In our simulation, the insulating TG-MHD to- 
tal energy spectrum can be well fitted globally only up 
to t = 2.2. After this time the energy spectrum displays 
a complicated behavior (see Fig. [4J. 

To study if this is an effect associated with insuffi- 
cient spatial resolution, in Fig. [5j we show the kinetic 
and magnetic energy spectra for the run performed on 
6144 3 points, as well as for a run with the same ini- 
tial conditions computed on a grid of 2048 3 points with- 
out bootstrap re-gridding, and as analyzed in [33]; we 
use both lin-log and log-log scales, for different times: 
t = 1.975, 2.201, 2.425 and 2.651. The implementation 
of the numerical procedure for the two runs in fact differs 
in several ways: (i) obviously, the resolution; (ii) single or 
double precision, the latter for the highest resolution; (hi) 
the truncation at high wavenumber (cubic for the latter, 



FIG. 5: (color online) Lin-log (a) and log-log (b) energy 
spectra as a function of wavenumber and for different times: 
t = 1.975, 2.201, 2.425 and 2.651, for the kinetic (light or 
red) and magnetic (dark or green) energy of the run at 6144 3 
resolution (the curves for this run go up to k max = 2048). 
Superimposed are the same for the ideal run in [33], at a res- 
olution of 2048 3 points in single precision (dark or blue for 
kinetic, light or brown for magnetic; the curves for this run 
go up to k m ax = 682). The same dominance of magnetic 
energy at small scales is observed in both runs; this implies 
extremely strong currents, compared to the vorticity, at small 
scales, as also observed when examining the temporal behav- 
ior of extrema (see Table [TT|. 



spherical for the former); and (iv) bootstrap regridding 
performed for the former, progressively in time. Yet, the 
two runs are seen to be equivalent. As time progresses 
in these flows, the magnetic energy gains from its kinetic 
counterpart (remember that Ey(t = 0) = E^it = 0)), 
particularly so at high wave numbers, as is also shown 
in Fig. [6] which gives the variation with wavenumber of 
the ratio EM(k)/Ev(k) for three different times and for 
both the 2048 3 and the 6144 3 runs. At t « 2.48, there is a 
surge of magnetic energy at small scales (large wavenum- 
bers) compared to its kinetic counterpart, a surge which 
finally resolves itself at the final time of the computa- 
tion. This behavior is likely linked to the evolution of 
structures in physical space (see { III C ) . 

When investigating the temporal evolution of the vor- 
ticity and current maxima, as shown in Fig. [7] (and also 
given in Table [B), we observe that there is a sudden 
change in the slopes at t « 2.5, and again at t « 2.65, 
the latter clearly discernible in the current density. These 
changes are associated with a shift of the maximum from 
one structure to another one. The first phase of evolu- 
tion, up to t « 2.5 is clearly exponential for both the max- 
ima of current and the vorticity, followed by faster growth 
on new structures that appear at later times. However, 
because the strongest peaks in current and vorticity ap- 
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FIG. 6: (color online) Ratio of magnetic to kinetic energy 
spectra for the I flow at three different times: t=2.49 (a), 
t=2.55 (b) and t=2.69 (c) for the run on a grid of 2048 3 points 
(circles, red) [33], and on a grid of 6144 3 (crosses, blue). The 
strong burst of excess magnetic energy at large wavenumbers 
subsides at later times. 



pear on a different structure at quite a late time in this 
run, when the grid resolution is almost reached, it is dif- 
ficult to ascertain whether a singularity would happen or 
not in this flow if it were pursued to yet higher resolutions 
and thus longer times. In other words, due to the physical 
structures that develop in this flow, the traditional tests 
of singularity (BKM and logarithmic decrement) cannot 
be applied in the latest evolutionary phase because it is 
too short. From that point of view, computations on yet 
higher-resolution grids will be necessary. In the next sub- 
section, we present a new analytical method that allows 
us to assess the plausibility of singularity scenarios. 



B. The link between the two known criteria for 
singularity 

It is known that several diagnostics for singularity can 
be used, and in fact that they are linked. The first 
method is to follow the temporal evolution of the maxi- 
mum of both vorticity and current and apply the BKM 
criterion given by Eq. M for fluids and generalized to 
the MHD case [21] (see [56] for the two-dimensional case 
in MHD); for smoothness on the [0, T] temporal interval, 
one must have convergence of the following integral: 



Using the Elsasser variables z ± = v ± b and defining the 
associated vorticities, U3 ± = u>±j, the above relation can 
also be written in characteristic form: 



"(.,f)||oo + ||w (.,t)|U) dt < oo 



Consider the formulation (111 for the BKM condition 



for regularity. If the numerical solution for the fields leads 
to a power-law behavior of the integrand, of the form 
1 1 w ( • , t ) 1 1 oo + 1 1 j ( • , t) 1 1 oo « C [T, - 1] -P , then the exponent 
/3 must be greater than or equal to one in order to be 
consistent with the existence of a singularity at time T* . 

The second tool for singularity diagnostic is to follow 
the logarithmic decrement d(t) of the fields mentioned 
above, in the context of the analyticity-strip method. In 
particular, one can look at the total energy spectrum (i.e., 
the spectrum of the sum of kinetic and magnetic energies) 
and calculate the decrement 6(t) for this spectrum. The 
logarithmic decrement 5(t) should go to zero in a finite 
time in order to be consistent with the existence of a 
singularity of the fields at time T*. In contrast, if S(t) 
decays exponentially in time then there is no evidence for 
a finite-time singularity. Finally, a third method consists 
of monitoring the evolution of the total production of 
small scales, through the enstrophy (integrated square 
vorticity) and the integrated square current. 

It may appear a bit odd to have different criteria to 
determine the evolution or not towards a singularity, but 
this is not redundant; quite the contrary. The link, at 
the level of heuristics, between the enstrophy divergence 
and that of vorticity was shown in pj [19] . More recently, 
a rigorous proof that bridges the two other criteria for 
singularity (BKM theorem and analyticity strip method) 
was shown in |20| along with an application to a numer- 
ical simulation of a 3D Euler fluid. The advantage of 
this bridge is that it leads to a new criterion when mon- 
itoring of the temporal evolution of small scales, giving 
an inequality between the power-law index of the energy 
spectrum and the temporal index of evolution for the 
logarithmic decrement, provided they can be assessed re- 
liably. 

To this end, one needs to use known inequalities. For 
our purposes, we recall the result in [20] that links the 
maximum vorticity modulus with the 3D Euler energy 
spectrum: 



(.,t)\\oo<cYk 2 y/E(kJ) , Vie[0,T) , (12) 



fc=i 



l«(..*)IU + ||j(.,*)||oo) dt< 



(11) 



where c is a constant of 0(1). 

The key concept in this new bridge is a hypothetical 
bound for the energy spectrum of the form 

E{k ) t)<Mk- nn{t) e- 2kSo{t) , Vie[0,T), V/seN, 

(13) 

for certain positive functions no(t) and 5o(t), and some 
positive constant M. The functions no,5o are closely re- 
lated to the analyticity-strip fit parameters Tlx , 5x con- 
sidered above, but they are not the same. In fact, the 
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above hypothetical bound is global (in fc-space), whereas 
as already mentioned the logarithmic decrement Sx(t) 
gives information on the asymptotic (large-fc) behavior 
of the energy spectrum. 

It was demonstrated in |20j that combining this hypo- 
thetical bound with the rigorous inequality ( 12 ) leads to 



a relation between the BKM theorem and the analyticity- 
strip method. To simplify matters, one considers the con- 
sequences of the following finite-time singularity scenario: 
suppose for simplicity that the exponent no in the hypo- 
thetical bound ( 13 ) remains constant as t approaches the 
singularity time T*, and that So(t) oc (T* — t) 1 , where 
7 > 0. Then the following necessary condition is found: 



1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.f 
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FIG. 7: (color online) Maxima of vorticity (brown dashed 
line) and of current (blue solid line) as a function of time for 
the I flow at high resolution. Note the jumps of the slopes 
(see also Fig. [9] below) corresponding to the emergence of 
different leading structures. 



7> 



6 - n 



in order that the blow-up be consistent with the BKM 
theorem. The formal argument is given in |20) and is 
immediately generalizable to MHD. The result for MHD 
is as follows: 



l|w(.,*)ll<»+l|j(. s *)IU < c 5>V2£ T (M) , Vte[o,T), 

fe=l 

(14) 

where now E T (k, t) represents the total energy spectrum, 
i.e., the sum of kinetic and magnetic energy spectra. The 
corresponding hypothesis for energy bound ( 13 ) is un- 



changed and similarly the hypothesis of blow-up for d T (t). 
The result is again a necessary condition, of the form 



7t > 



G 



(15) 



Note that, since in the Euler case, the observed uq ap- 
pears to be (at least for some initial conditions) larger 
than the exponent n T in the MHD case, one sees that 
the eventual realization of a singularity in MHD might 
be a different process than for the Euler equation. This 
is not necessarily surprising for at least three reasons: (i) 
MHD is thought to be smoother than hydrodynamics, 
insofar as Alfven waves may slow down the dynamics of 
propagation to small scales, leading possibly to a different 
energy spectrum, the so-called Iroshnikov-Kraichnan law; 
(ii) the Onsager principle concerning energy dissipation 
can likely be replaced in MHD by magnetic helicity con- 
servation, following the so-called Taylor conjecture [21], 
thereby changing the dimensionality of the system; and 
(hi) the degree of smoothness required to ensure total 
energy conservation (technically, the index of the Besov 
space needed) for the velocity and the magnetic field may 
differ in a way that is compatible with the Iroshnikov and 
Kraichnan spectra |21j . In particular, with uq rs 4 for Eu- 
ler, one obtains 7 > 1 whereas for n T s» 3 in ideal MHD 
(see [32]), one has 7 T > 2/3: the decay of the logarithmic 
decrement would be slower in MHD, as expected because 
of the slowing-down of the dynamics by (Alfven) waves. 



1. Analysis of the total energy spectrum 

At time t w 2.33, we observe a change in the behavior 
of the total energy spectrum, probably due to the im- 
minent, accelerated collision between two current sheets 
(confirmed by inspection of the structures in real space) , 
and the corresponding fast generation of a second length 
scale, related to the distance between the two sheets. The 
original length scale of the problem, interpreted as the 
decreasing width of the current sheets, decreases slower 
than this new length scale so eventually the two length 
scales become comparable. It is known that when two 
or more sharp physical structures of similar length scales 
are present, the traditional fit (10) of the energy spec- 



trum fails. For example, in the Kida-Pelz 3D Euler flow, 
the departure of the measured energy spectrum from the 
traditional form ( 10 1 was modeled with good accuracy by 



attributing it to interferences of two complex singularities 
situated at equal distances from the real axis [54 . How- 
ever, the extra complexity (spatial and temporal) of the 
MHD flow under current study makes it difficult for us to 
find a good model for this new behavior. This imposes a 
practical limitation on the analyticity-strip method as a 
means for finding a good estimate of the actual logarith- 
mic decrement <5 T (t) of the spectrum (where "actual" is 
used in contrast to the measured one). In fact, depend- 
ing on the fit interval we get vastly different estimates 
for the width 6 T (t) for times t > 2.33, so our knowledge 
of the width 5 T (t) as in the large-fc asymptotic expansion 
\nE T (k,t) ~ — 2kS T (t) has significant errors that grow 
in time. 

In conclusion, we cannot tell by using the analyticity- 
strip method alone whether there is a finite-time 
singularity in the MHD flow under study at times 
t > 2.33. Of course, we know from continuity arguments 
that the width 8(t) should remain non-zero at least for 
a short time after t = 2.33. But that is all we know, so 
there are two possible scenarios: 



Scenario 1: There is no finite-time singularity up to 
time t = 2.7, so the simulation is well resolved, 
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FIG. 8: (color online) Running estimates of the analyticity- 
strip method exponents n T (t) (a) and S T (t) (b) for the to- 
tal energy spectrum. In dashed (blue) line for a fit interval 
[4, 2048] ; in thick solid (magenta) line for a fit interval [4, 500] ; 
in thin with dots (brown) line for a fit interval [500,2048]. 
In (b), the horizontal line represents the reliability threshold 

^T^max = 2. 
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FIG. 9: (color online) a) Sum of maxima of vorticity and 
current as a function of time. There is a clear jump at t = 
2.476 corresponding to an emergent near-singular structure 
taking over a previous one. b) Multiplicative inverse of the 
logarithmic derivative of the previous curve. If this has a 
negative slope, it is an indication for a possible finite-time 
singularity. 



perhaps marginally. The implications of Scenario 
1 will be exploited in Section fill C| 



Scenario 2: There is a finite-time singularity at a time 
between t = 2.33 and t = 2.7, but this cannot be 
assessed using the analyticity-strip method alone. 



Let us consider the implications of Scenario 2. Although 
we do not know the logarithmic decrement r5 T (£), we can 
still have an estimate for the positive exponent n T (t) ap- 
pearing in the bound ( |13[ ) for the total energy spectrum. 
In fact, what is needed in inequality ( |15[ ) is a lower bound 
for n T (t), rather than n T (t) itself. This lower bound can 
be estimated by looking at the low wavenumber fits of 
the total energy spectrum (from k = 4 to k = 500), as 
shown in Fig. [8] Running estimates for n T (t) obtained 
in this way turn out to be consistently smaller than the 
estimates obtained by using fit intervals including larger 
values of k. The result for the lower bound is n_ = 2.385. 
With this number, the inequality ( |15[ ) gives a bound for 
the unknown exponent 7 T in S T (t) oc (T* — t) 7T : 

7t > tt^— « 0.553 , (16) 
6 — n_ 

so even though we do not know whether the logarithmic 
decrement is going to zero or not in a finite time, we 
have been able to estimate how fast it should go to zero 
in the hypothetical case of a finite-time singularity. 



2. Analysis of the sum of supremum norms of vorticity and 

current 

To further comment on the feasibility of Scenario 2, 
let us consider the method of running estimates for sin- 
gularity of fast-growing quantities introduced in |19j . 
We apply this method to the growth of the BKM field 
I !<*>(., t)| |oo + | *) | loo with the ansatz ||u>(., *)!!«, + 
||j(-,*)||oo ~ C[T* -t]~' 3 . The method gives running es- 
timates of the exponent (3 and of the singular time T*. 
In Fig. [9]ja) we observe that there is a jump at t = 2.476 
in the growth rate of the BKM quantity; however, this 
is not due to a dynamical effect. It is rather due to an 
independently emergent physical structure that is more 
singular than the previous one. Figure [9jb) shows the 
multiplicative inverse of the logarithmic derivative of the 
BKM quantity. If this curve has negative slope, then 
the intersection of the slope with the i-axis gives a run- 
ning estimate of the potential singularity time. We see 
two instances of negative slope. We discard the instance 
at about t = 2.476 because this is due to the transient 
emergence of the new structure. However, near t = 2.5 
we observe more robust evidence of potential singularity, 
although the data is quite noisy and thus the tangent is 
oscillating too much, so we cannot have precise estimates 
of the singularity time and the exponent /?. Naked-eye 
prediction of singularity time, obtained by finding the in- 
tersection of the smoothed tangent at t — 2.5 with the 
t-axis, would give T* sa 2.56-2.58 and /3 « 1.44-1.75. 
These values for the estimates of T* and /3 were obtained 
by estimating two tangents in Fig.[9](b), each tangent be- 
ing defined as the linear interpolation of six contiguous 
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TABLE II: Time, maxima of current, and their kj) lo- 

cation in the fundamental [0, 7r/2] box in grid units, as well as 
maxima of vorticity and their location, for the high-resolution 
I flow on a grid of 6144 3 points. The indices (i, j, k) refer to the 
grid points in (x, y, z) where the maxima take place. Note the 
sudden jumps in the position of the maxima; the first jump in 
coordinates and in values of maxima occur for t ~ 2.48, and 
the second one at t ~ 2.62. 



data points taken out of the seven data points highlighted 
in the figure. 

It is interesting that near t = 2.5 the estimated loga- 
rithmic decrement S T (t) (using the full fit range [4, 2048]) 
indeed has a change in behavior, first a deceleration and 
then an acceleration, although this occurs near the reli- 
ability threshold-see Fig. ^b). A computation of the 
running estimate of decay exponent 7 T as in S T (t) oc 
(T* -t) 7T gives 7 T = 0.94 at t = 2.501 but only that data 



point agrees with the rigorous bound in ( 16 ), 7 T > 0.553. 
At slightly later times, the estimated value of j T becomes 
10 times smaller, thus violating the rigorous inequality. 
The corresponding predicted singular time, using this 
method, gives a running estimate « 2.516-2.522. 

To summarize, Scenario 2 is plausible but some of its 
aspects occur in the limit of the reliability threshold. 
This point is aggravated by the fact that the sampling of 
current and vorticity maxima at the grid points induces 
spurious oscillations in the data (a way to suppress these 
oscillations is discussed in Sec 
conclusion can be drawn at 
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the 



Therefore, 
moment. 



no robust 
A future 



higher-resolution numerical simulation should shed more 
light on the feasibility of Scenario 2. 



C. Structures in physical space 

Visualization plays an important role in the discovery 
process, and many of the arguments considered above 
used information from the evolution of the structures 
in physical space, based on previous runs |32j and the 
present high-resolution computation. In order to visual- 
ize the velocity and magnetic field and their gradients, 
one needs to reconstruct the three-dimensional data us- 
ing the four-fold symmetries of the TG-MHD configura- 
tion, a daunting task at such resolutions. In that con- 
text, note that the VAPOR visualization system devel- 
oped at NCAR j55[ [56] allows one to analyze the data 
using wavelet compression in order to explore rapidly at 
coarser resolutions, and then to increase the resolution 



as needed where needed. 

The acceleration in the formation of small scales was 
first identified in |32j with the collision of two current 
sheets leading to a quasi-rotational discontinuity. The 
present computations at higher resolution confirm these 
results and allow us to go further in time. We have given 
in Table [Tl] the values close to the end of the compu- 
tation of the maximum of the vorticity and of the cur- 
rent as well as their location in the fundamental [0,7r/2] 
computational box. Concentrating on the current, which 
is known in two dimensions to have a simpler geometric 
structure (a dipole instead of a quadrupole for the vortic- 
ity), we observe two jumps, near t\ = 2.48 and ti = 2.62, 
both in the value of the current maximum and in its lo- 
cation. 

The collision of two current sheets leading to a quasi- 
rotational discontinuity was clearly observed in [32] at 
a resolution of 2048 3 points; this phenomenon is con- 
firmed in the present computation with three times the 
linear resolution and thus any numerical effect can be 
ruled out for it. The second acceleration in the develop- 
ment of small scales, which occurs at a later time, seems 
to be related to the near co-location of these two sheets 
and this is now what we examine by considering struc- 
tures in physical space. We also note that such current 
sheets are known to roll-up at sufficiently high resolution 
in the dissipative case, and similar rolled-up structures 
have been observed in the Solar Wind in a much more 
complex physical environment [571 158]. But they are only 
a recent finding in DNS of MHD turbulence on grids of 
1536 3 points using the GHOST code [MSI] and 2048 3 
(equivalent) points using TYGRS [4"Tj . 

The maximum of the current first moves along the ver- 
tical axis (index k), then at t = 2.5 it has moved in the y 
direction (index j): it is traveling along the lower sheet 
as the two sheets seem to join with each other, to follow 
the curvature. Then finally, around t — 2.65, the current 
density maximum now moves along the diagonal in the 
horizontal plane. 

A rendering of the current is given in Fig. [I0]at t=2.54 
(left and middle for the vorticity and the current). It is 
the merging of two current sheets that causes the maxi- 
mum to go radially (in a cylindrical sense) from the cor- 
ner of the fundamental box ([i, j, k] = [1, 1, 1537]) along 
a polar angle of 7r/2 on the top plane. The two sheets 
(seen at both ends of the box because of symmetries) 
are clearly almost touching each other (a zoom indicates 
that they are two to three grid-points apart). Finally, on 
the right is given a suite of six magnetic field lines at the 
latest time of the computation that appear to all con- 
verge to one point, indicative of a potential singularity, 
that point being the location of the current maximum 
at that time. The two current sheets are barely visible 
(purple, and blue below). Alternative views are given in 
Fig. [Tl] with in particular two-dimensional cuts at the 
highest (grid) resolution, indicating that the two current 
and vorticity sheets are still individually resolved. 

Such features correspond to a strong bending of mag- 
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FIG. 10: (color online) Perspective volume rendering using the VAPOR software [55, 56 of the vorticity (left) and of the 
current density (middle) at t=2.54. Note the occurrence of a double layer structure due to the collision and subsequent joining 
of two sheets. At a later time (t=2.65; right), the magnetic field lines taken on these two colliding sheets all go to the same 
location, which coincides with the maximum of the current (coordinates given in the Table), implying sharp localized bending 
(and possibly torsion) of magnetic field lines in the vicinity of that maximum. 




FIG. 11: (color online) Alternative views of the same struc- 
tures at the same time as in Fig. |10| two-dimensional cuts 
of vorticity (left) and of current (right) in the region of the 
current sheet collision, displayed down to the grid-resolution. 



netic field lines in the vicinity of the current and vorticity 
maxima, implying strong directional variations. It may 
also imply magnetic field line stretching in this strong 
curvature region, a stretching that would be consistent 
with the sudden increase in magnetic energy (relative to 
its kinetic counterpart), as observed clearly (see Fig. [6]). 
This is also reminiscent of the necessity, in the Euler 
case, of a blow-up of both the magnitude of the small- 
scale field but also of the curvature of its field lines, as 
shown in [62 , 63 , for a singularity to occur. 



D. The case of other Taylor-Green configurations 
in ideal MHD 



Finally, let us mention briefly how the two other initial 
conditions satisfying the TG symmetries behave, with 
ideal runs computed on grids of up to 4096 3 points. Sim- 
ilar temporal evolutions seem to occur for both flows, as 
shown in Figs. fl2] and |13) The spectral indices seem to 
reach values smaller than in the Euler (ideal fluid) case 
for all configurations examined in this paper, systemat- 
ically below a value « 3, with some oscillations in the 



conducting case (C flow, Fig. 12 two lower panels). 



The maxima of current and vorticity are displayed in 
Fig. [13] for the A and C configurations. The C flow cur- 
rent maximum (Fig. 13 center) undergoes first a jump 
from structure to structure at relatively early times, fol- 
lowed by a traditional exponential phase corresponding 
to the thinning of current. This is followed again by a 
short and rapid further increase in the maximum which 
appears difficult to analyze in more detail, due to the fact 
that the temporal interval during which this latest ac- 
celeration occurs is again too short, as for the insulating 
configuration analyzed in the previous sections. The vor- 
ticity maximum in the A flow shows a rather monotonic 
increase until the final acceleration in current. For the A 



flow (Fig. 13 top) the monotonicity of the current and 
vorticity maxima are reversed compared to the C flow. 
We note that the temporal evolution of the integrated 
square vorticity and current for the A and C runs indicate 
that they become nearly equal for late times (not shown) . 
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This is due to the fact that, after the grid resolution is 
reached by the velocity and magnetic field structures, the 
evolution is that of a truncated system of Fourier modes 
which evolve, in the simplest case, to equipartition due 
to statistical equilibrium, as analyzed in [64] ; this begins 
to occur at late times in these computations, at a faster 
rate the smaller the scale. 

Another instance of quasi-equipartition between ki- 
netic and magnetic energy is occurring at earlier times, 
and is reminiscent of what is observed in the dissipa- 
tive case for many configurations (see, e.g., [33j|42]). In- 
deed, we see that the ratio of the spectra of magnetic 
and kinetic energy given in Fig. [13] (bottom) for the C 
configuration is close to (and slightly above) unity from 
k ~ 4 up to the maximum wavenumber. This is also 
observed for the other two configurations examined in 
this paper and is consistent with the expression for the 
spectra obtained for ideal dynamics of a truncated sys- 
tem with zero (or negligible) helicity [64] . in which case 
equipartition obtains. This appears to be another exam- 
ple where the ideal dynamics is consistent (and can be 
viewed as predictive of) dissipative (and/or forced) iner- 
tial range dynamics, as first clearly showed using direct 
numerical simulations in the fluid case in [3D]. We also 
note that such a quasi-equipartition of kinetic and mag- 
netic energy, with in most cases a slight excess of the 
latter, has been observed for a long time in Solar Wind 
data |65j , and confirmed later by more detailed observa- 
tions as well. 



IV. SUMMARY OF RESULTS AND 
CONCLUSIONS 

We have shown in this paper several new results 
concerning the ideal dynamics of MHD configurations, 
namely that: (i) by increasing the resolution by a fac- 
tor of three from our previous study, we still reproduce 
the results obtained in a 2048 3 simulation up to the last 
time computed in that run, including an acceleration in 
the maximum of current and vorticity and in the de- 
crease of the logarithmic decrement; (ii) in the new high- 
resolution simulation, we see yet a second acceleration of 
the formation of small scales at a later time, in a situation 
that is as well resolved as the previous acceleration was 
in the 2048 3 simulation; (iii) these two accelerations are 
clearly associated with changes in the structures in phys- 
ical space of the current and vorticity; (iv) these changes 
also pollute the small-scale spectrum, creating a limita- 
tion in practice to the applicability of the analyticity strip 
method; (v) a new method, bridging the analyticity strip 
method and the so-called BKM criteria by means of sharp 
analysis inequalities, is extended to MHD, and allows us 
to rule out spurious singularities; (vi) the new method 
cannot completely rule out the existence of a finite-time 
singularity at a time between t = 2.33 and t = 2.7; (vii) 
the structures that seem to create this acceleration in 
the formation of small scales are related to the near colli- 
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FIG. 12: (color online) Analysis of the fit to the total energy 
spectra for the "A" flow (two upper panels) and "C" flow (two 
lower panels) on grids of up to 4096 3 points; note that early 
times are not shown. First and third panels: Logarithmic 
decrement; the dashed line indicates the resolution limit for 
the 4096 3 grid. The smallest grid resolution is reached at 
t « 1.85 for the A flow and at t ~ 2.35 for the C flow. Second 
and fourth panels: Spectral index of the total energy, with 
values that seem to settle below n — 3 for both configurations. 



sion and further spatial co-location of two current sheets, 
and similarly for the vorticity; (viii) these results do not 
seem to be occurring only for one flow, but seem to take 
place as well in the other two configurations studied in 
this paper, up to equivalent resolutions of 4096 3 points; 
and finally (ix) a simple re-gridding technique, which al- 
lows for substantial savings in computer time, is shown 
to be entirely reliable provided a conservative threshold 
for applying the method is utilized. We should note that 
in one case (that of the I configuration) , the scale sepa- 
ration reached in the computation is unprecedented up 
to this point in time. 

In summary, we have found that at high resolution, the 
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most intense structures that develop in ideal MHD come 
from the near collision and later from the near juxtapo- 
sition of two current and vorticity sheets. The maxima 
of these small-scale fields undergo abrupt jumps twice, 
and it will be necessary to pursue this computation at 
yet higher resolutions to see whether the criteria for a 
singularity to develop or not are satisfied, by monitoring 
for a time that is sufficiently long the maxima of current 
and vorticity and to compute other diagnostics as well. 

We remark that our pursuit is not just a brute- force in- 
crease in resolution. In fact, we have made use of a new 
analytical tool, that bridges two known singularity cri- 
teria (BKM-type theorem for MHD and analyticity-strip 
method) , leading to a new method for ruling out spurious 
indications of singularity. We have applied this method 
to the current configuration under study at the highest 
resolution achieved in this paper, and concluded that the 
existence of a finite-time singularity at a time between 
t = 2.33 and t — 2.7 cannot be completely ruled out. 
While it would be desirable to produce a more specific 
statement in this regard, there is one fact that makes 
it difficult to advance further: the values of £)||oo 
and ||j(-, <)||oo, needed for testing singular behavior in the 
framework of the BKM theorem are currently measured 
using collocation-point data, a standard procedure that 
leads to "noise" or error in the data. This noise is evident 
as tiny oscillations in Fig. ^b) which add an uncertainty 
to the computation of slopes. The source of this noise was 
discovered recently in |66) , in the context of the more con- 
trollable inviscid Burgers one-dimensional flow. There, 
as in our MHD case, the systematic periodic sampling of 
collocation-point maxima introduces an error in the pre- 
cision of the measurement with respect to the true value 
of maxima. The error oscillates in time; its frequency 
grows with the numerical resolution used if the time step 
is determined by a fixed-ratio CFL condition. Moreover, 
the amplitude of the error depends on the spatial profile 
of the maximum computed, so that the error increases 
as the structures become more peaked. In [66] the solu- 
tion to this problem was proposed and has two levels of 
complexity: at the simplest level, a post-processing com- 
putation of extrapolated values of the maxima of vor- 
ticity and current can eliminate partially the oscillatory 
part of the error. At the deepest level, the application 
of an adaptive time stepping beyond CFL, so that the 
product At x (||u;(-,t)|| 00 + ||j(-, t)||oo) remains constant, 
can improve the precision in the computation of vorticity 
and current maxima by a factor of 10 2 , at no extra cost in 
computational time and memory [(JB] . In our future work 
we will implement these procedures so we can have more 
robust evidence regarding the hypothesis of finite-time 
singularity in MHD. 

One dynamical effect that can play a role in stopping a 
putative singularity is the phenomenon of dynamic align- 
ment that is rather ubiquitous in turbulent flows. For 
example, it was shown in [43j that the alignment of vor- 
ticity with shear or pressure gradients, and equivalently 
of magnetic field and shear, enhances point- wise helicity 
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FIG. 13: (color online) Top, Center: Maxima of current 
(dashed, red) and vorticity (solid, black) over time for the 
A configuration (top), and for the C configuration (center). 
The abrupt changes at intermediate times are linked to the 
fact that the maxima jump from structure to structure. Note 
the rapid increase at t > 1.85 in the C flow and, and t > 2.35 
for the A flow, and the noise at late times due to the fact that 
small-scales have become insufficiently resolved. Bottom: Ra- 
tio of the magnetic to kinetic energy spectra, for the C config- 
uration at different times, with dark blue (0.43 < t < 0.79), 
green up to t — 1.15, red up to t — 1.51, black up to t — 1.87, 
magenta up to t — 2.23, and cyan up to t = 2.48. Note the 
quasi-equipartition in the inertial range all the way to the 
cut-off, as well as the excess magnetic energy at large scales. 



15 



(kinetic helicity in the former case, cross helicity in the 
second case), although the global norms are conserved, 
and it does so in a time of the order of the eddy turn- 
over time. In fact, an alignment between all variables 
involved in the nonlinear terms of MHD, namely, veloc- 
ity and magnetic field in Ohm's law [67], velocity and 
vorticity in the Lamb vector, and current and magnetic 
field in the Lorentz force, occurs rather systematically, in 
particular the latter [68J . It is not clear what the effect of 
dissipation is in these alignment properties, or whether 
such alignment tendencies would be sufficient to prevent 
singularities to occur in the ideal case. In that light, a 
more detailed analysis of the local properties of the flow 
in the vicinity of the current and vorticity maxima will be 
undertaken in a follow-up paper. Furthermore, ideal and 
dissipative flows have common properties because of their 
nonlinear multi-scale interactions. The lack of universal- 
ity, found in decaying flows with imposed Taylor-Green 
symmetries [33J is also found in the forced case [59], and 
thus it is an open problem to see whether it will occur in 
the ideal case, although the differences between inertial 
indices is small and thus requires high resolutions and 
long-time integration. 

A theory of turbulent flows is still lacking, and yet 
such flows are ubiquitous in nature and are an integral 
part of the problem of weather prediction, of climate as- 
sessment, of understanding the formation and prediction 
of extreme events such as tornadoes and hurricanes, of 
reconnection events in space physics such as solar flares 
and coronal mass ejections, plasmoids, and in disruptive 
plasmas. Such flows develop intense small scale struc- 
tures in the form of vortex and current sheets and fila- 
ments with power-law scaling properties and departure 
from Gaussianity attributed to intermittency. Similarly 
in the ideal case at intermediate times and intermediate 



scales, a classical turbulent spectrum has been observed 
recently for fluids [401 ED] , with at smaller scales the sta- 
tistical equilibrium that can be derived analytically using 
the quadratic invariants preserved by the truncation (see 
[64] for 3D MHD), the whole flow evolving progressively 
towards flux-less Gaussian equilibrium solutions. What 
is lacking is, among other things, a statistical description 
of the small scales, and a prediction of long-time large- 
scale dynamics with ensuing modified transport proper- 
ties. By combining this study with a well-resolved high 
Reynolds number dissipative run, one may be able to es- 
tablish in 3D-MHD the link between the role of ideal 
nonlinear dynamics, and dissipative- induced reconnec- 
tion (see e.g. [H]), leading to finite dissipation in the 
limit of zero viscosity and magnetic resistivity as shown 
in both two-dimensional [711 172j and three-dimensional 
cases [61] . This may shed light on dissipation processes in 
turbulent conducting flows, and on the role of non-local 
interactions between disparate scales [73J [73] in MHD 
when compared to the Euler (fluid) case (see also [75]). 
thus leading to better estimations of the energy dissi- 
pation rate controlled by turbulence in astrophysics and 
space physics. 
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